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Summary. Taking into account that a proper description of disordered systems 
should focus on distribution functions, the authors develop a powerful numerical 
scheme for the determination of the probability distribution of the local density of 
states (LDOS), which is based on a Chebyshev expansion with kernel polynomial 
refinement and allows the study of large finite clusters (up to 100 3 ). For the three- 
dimensional Anderson model it is demonstrated that the distribution of the LDOS 
shows a significant change at the disorder induced delocalisation-localisation transi- 
tion. Consequently, the so-called typical density of states, defined as the geometric 
mean of the LDOS, emerges as a natural order parameter. The calculation of the 
phase diagram of the Anderson model proves the efficiency and reliability of the 
proposed approach in comparison to other localisation criteria, which rely, e.g., on 
the decay of the wavefunction or the inverse participation number. 



1 Introduction 

The localisation of quantum particles in disordered systems is one of the most inten- 
sively studied problems in condensed matter physics [1-5]. In real systems disorder 
can arise for a number reasons. We may think of randomly distributed impurities, 
vacancies or dislocations in an otherwise ideal crystal, of random arrangements of 
electronic or nuclear spins, etc. While the disorder appears in many forms that are 
sometimes difficult to characterise theoretically, the randomness in the model intro- 
duced and discussed by Anderson is simple but sufficient to capture the basic features 
of the disorder- induced metal insulator transition [6]. The Anderson Hamiltonian, 



describes noninteracting electrons moving on a lattice with random on-site potentials 
(compositional disorder). The operators cj (cj) create (annihilate) an electron in 
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a Wannier state centred at site j, and the local potentials tj are assumed to be 
independent, uniformly distributed random variables, 



The parameter W is a measure for the strength of disorder and is usually given 
in units of the nearest neighbour hopping matrix element t. Throughout this work 
we consider d-dimensional hyper-cubic lattices with N = L d sites, impose periodic 
boundary conditions (PBC), and set the lattice spacing equal to unity. 

The spectral properties of the Anderson model (1) have been carefully analysed 
(see, e.g., Ref. [7]). For sufficiently large disorder or near the band tails, the spectrum 
consists exclusively of discrete eigenvalues, and the corresponding eigenfunctions are 
exponentially localised. Since localised electrons do not contribute to the transport 
of charge or energy, the energy that separates localised and extended eigenstates is 
called the mobility edge. For any finite disorder W > 0, on a one-dimensional (ld) 
lattice, all eigenstates of (1) are localised [8,9]. This is believed to hold also in 2d, 
where the existence of a transition from localised to delocalised states at finite W 
would contradict the one parameter scaling theory [10,11]. 

In spite of the progress made over the last decades, the Anderson metal-insulator 
transition is still not completely understood. There are several reasons why the ex- 
isting theories remain unsatisfactory. Especially when electron-electron or electron- 
phonon interactions come into play, the very successful one-parameter scaling ap- 
proach might be problematic, because close to the localisation transition the energy 
scales associated with both disorder and interactions are comparable to the Fermi 
energy [12]. On the other hand, the numerical study of the localisation-delocalisation 
transition is demanding, since the involved length scales can become extraordinary 
large, in particular near the critical point. Obviously, methods that are based on a 
full diagonalisation of the Hamiltonian and on the study of the one-particle eigen- 
states are restricted to rather small systems. Examples are the calculation of the 
localisation length from the decay of the electronic wavefunction or the evaluation 
of the inverse participation number. In addition, one-particle eigenstates are not de- 
fined for interacting systems. Hence, none of these criteria can easily be generalised 
to interacting disordered systems. To overcome these difficulties is perhaps the most 
challenging issue of current research on disordered materials. 

Motivated by this situation, this contribution provides a (quasi approximation 
free) numerical analysis of the recently revived "local order parameter" approach to 
the Anderson transition, which, within the framework of the statistical dynamical 
mean field approximation, has been successfully applied also to correlated electron 
(phonon) systems [12,14]. Adopting a local point of view and focusing on the distri- 
bution of the physically interesting quantities, the method follows the original route 
to the localisation problem established by Abou-Chacra et al. [16]. In particular, 
we demonstrate that for Anderson type models the distribution of the local density 
of states (LDOS) can be determined very easily by the kernel polynomial method 
(KPM) [17], a refined Chebyshev expansion technique. Based on the distribution of 
the LDOS, localised states are distinguished from extended states by a vanishing 
geometrical average, which is usually called the "typical DOS". In addition, it turns 
out that this quantity characterises the disorder-induced metal-insulator transition 
also in more complex systems [12, 14, 19-22]. 
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To examine the efficiency and accuracy of the proposed LDOS-KPM approach, 
we carry out a comparative numerical study of the localisation-delocalisation tran- 
sition. In view of the wealth of known results, the 3d Anderson model seems to 
be best suited for this kind of investigation. The results we obtain for the mobility 
edge from different methods allow for a detailed understanding of the typical DOS 
concept and open the road towards an application to more complex situations. 



2 Anderson transition 

As the Anderson transition is expected only for d > 2, in this chapter we focus 
on the 3d case, for which the lack of successful analytical approaches necessitates 
a numerical treatment. In contrast to the widely used numerical transfer matrix 
method [23-25], which describes the 3d system as a quasi-ld system of finite cross 
section, below we stress the bulk properties of the system and consider cubic clusters 
which extend equally in all spatial directions. Due to the large length scales that 
emerge in the critical region, it is generally a difficult task to interpret the results 
of such finite cluster calculations. 

As a kind of benchmarking, we review and compare established localisation cri- 
teria, namely, the localisation length (Sec. 2.1) and the inverse participation num- 
ber (Sec. 2.2). Both can be extracted from the one-particle wavef unctions, which, 
however, requires the complete numerical diagonalisation of the Hamiltonian (1). 
In Sec. 2.3, we present the new approach that is based on the distribution of the 
LDOS. Since the calculation of the LDOS via the KPM requires only sparse matrix 
vector multiplications, this technique scales linearly with the number of lattice sites 
and permits the study of significantly larger systems. 



2.1 Decay of the wavefunction 



W=3t 
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Fig.l Decay of an elec- 
tronic wavefunction ip n in 
the band centre as a func- 
tion of the distance r = 
|r< — r max | to the site 
with maximum amplitude, 
r m ax. Results are given 
each for one fixed energy 
and realization of disorder, 
N = 30 3 , PBC. 



The most obvious but costly way to access the localisation properties of an 
electronic wavefunction is the direct calculation of the localisation length A, which 
is infinite for extended states and finite otherwise. For localised states, the envelope 
of the wavefunction decays exponentially from some point r max in the crystal. 
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|Vn(Pi)| ~ /(rt)exp A rmax ' j . (3) 

The random function /(rj) describes the statistical fluctuations of the amplitudes 
V'n(ri) of the eigenfunction ip n at energy E„. Given ij) n , the localisation length X(E n ) 
is obtained by locating the site of maximum amplitude, r max , and fitting Eq. (3) to 
the data. In contrast to the case of weak disorder, where the amplitude is essentially 
independent of the distance from r max , at higher values of W a clear exponential 
decay is observed (see Fig. 1). Note, that besides the direct fit with equal weight for 
the amplitudes of all sites, A can also be determined using the method of asymptotic 
slope [26] . Here the data is first averaged within shells of fixed distance from r max and 
fitted thereafter. However, using this second approach the detection of the Anderson 
transition is not as robust and more sensitive to the fluctuations of the data. We 
therefore refrain from considering corresponding results. 



2.2 Inverse participation number 




Fig. 2 Upper part: Averaged 
inverse participation number 
Pav 1 of 1000 systems with L s 
sites and PBC (top to bottom: 
L = 8,9,10,12,14,16,20). 
Lower part: Probability den- 
sity of P" 1 for the 10 3 sys- 
tem in the band centre (solid 
line) and near the band edges 
(dashed line). Note the differ- 
ent scales in the lower panels. 



Yet another quantity that measures the Anderson transition is the inverse par- 
ticipation number [27], 

N 

P-\E U ) = Y,\i>n(Ti)\\ (4) 

i=l 

which is proportional to the inverse number of sites that contribute to a given one- 
particle wavefunction ip n . For delocalised states we find P _1 ~ 1/JV, which vanishes 
in the thermodynamic limit. Localised states, on the other hand, approximately 
extend over a finite volume iVo, yielding P~ ~ 1/JVq independently of the system 
size N. In Fig. 2 this behaviour is demonstrated for small and large disorder W. 
While in the localised case (W = ISt) P" 1 is almost independent of L, apparently 
it decreases with L for extended states (W = Wt). Apart from the different scaling, 
the distribution of P~ l changes at the transition, being sharply peaked for extended 
states and rather broad for localised ones (lower part of Fig. 2) . 

Based on the distribution of the participation number P recently an alternative 
numerical approach for monitoring the Anderson transition was proposed [28]. In 
analogy to results for a certain class of power-law random banded matrices, which 
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Fig. 3 Normalised standard 
deviation of the participa- 
tion number in the band cen- 
tre as a function of disor- 
der for different system sizes 
using PBC. The obtained 
results were averaged over 
1000 realizations of disorder. 
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indicate the scale invariance of the distribution of P at the Anderson transition, 
Malyshev et al. [28] suggest to detect the transition by studying the ratio of the 
standard deviation of P, AP, to the mean participation number P av , which should 
be independent of the system size at the critical disorder. 

So far this approach has only been tested for a one-dimensional model with 
diagonal disorder and power-law long-range hopping [28], which shows a transition 
at the band edge and can thus be tackled with the Lanczos method. In Fig. 3 
we present first data for the band centre of the standard Anderson model. While 
for small disorder the ratio AP/P^ decreases with increasing system size, at large 
disorder the opposite happens. The intersection of the curves is not completely 
independent of the system size (see inset of Fig. 3), and a precise determination of 
the transition requires some finite size scaling of the data. Performing a finite-size 
scaling [28] our data is consistent with a critical disorder strength of W c « 16.1 ±0.8 
in the thermodynamic limit. 

2.3 Local density of states 

Already in his seminal paper [6] Anderson pointed out that in order to describe 
the transition from delocalised to localised states it is very instructive to discuss 
the distribution of local quantities of interest, such as the escape rate or recurrence 
probability from or to a given site. Another suitable quantity that becomes critical 
at the Anderson transition is the LDOS [12,29], 



which for a given energy directly measures the local amplitude of the wavefunction 
at site Vi. So far the LDOS has been considered mainly within analytical approaches 
or by the Lanczos recursion method [29]. Typically the latter suffers from severe 
stability problems at high expansion order and conclusive results for the Anderson 
transition are difficult to obtain [30] . Fortunately the KPM technique [17] described 
in Appendix A is a very efficient way to circumvent these difficulties and allows the 
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Fig. 4 General shape 
and finite size scaling of 
the LDOS distribution 
p(pi/pav)- Keeping the 
ratio N/M = 1.95 fixed 
for N = 10 3 , 20 3 , 40 3 
and K r x Jf s = 
10 4 xl00, 100x100, 32x32 
respectively, we calcu- 
lated histograms for 
E e [-0.lt, 0.1*]. Inset: 
Double logarithmic plot of 
p(Pi/Pav) for the localised 
case together with a 
log-normal fit to the data. 



calculation of high-resolution LDOS data for very large systems. In a nutshell, within 
this approach the function of interest is expanded in a finite series of Chebyshev 
polynomials. To weaken the effects of the truncation and ensure properties such as 
positivity and normalisation, the function is convoluted with an appropriate integral 
kernel. The resolution of the method is inversely proportional to the order of the 
expansion M (the number of so-called Chebyshev moments). 

Adopting Anderson's original point of view that a proper description of disor- 
dered systems should focus on distribution functions, we calculated pi{E) for a large 
number of samples, K r , and sites, K s , and studied its statistical properties. In Fig- 
ure 4 we show the resulting distribution of pi(E = 0), normalised by its mean value 
p av , for two characteristic values of disorder. As p av is a function of disorder, this 
normalisation ensures (pi/p av ) = 1 independent of W, allowing thus an appropriate 
comparison. In the delocalised phase, W = 3i, the distribution is rather symmetric 
and peaked close to its mean value. Note that increasing the system size and the 
expansion order, such that the ratio of mean level spacing and KPM resolution is 
fixed, does not change the distribution. This is in strong contrast to the localised 
phase, e.g., W = 24i, where the distribution of pi(E) is extremely asymmetric. Al- 
though most of the weight is now concentrated close to zero, the distribution extends 
to very large values of pi, causing the mean value to be much larger than the most 
probable value. In addition, a similar finite size scaling increases the asymmetry and 
underlines the singular behaviour expected in the thermodynamic limit and at infi- 
nite resolution. Note also, that the distribution of the LDOS is well approximated 
by a log-normal distribution [31], 

as illustrated in the inset of Fig. 4. 

Of course, the study of entire distributions is a bit inconvenient, and for practical 
calculations, instead, we will prefer an appropriate statistics that uniquely charac- 
terises the distribution. The above findings suggest, that such a statistics is given 
by the arithmetic and geometric averages of pi{E), 
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Fig. 5 Average 
(solid line) and 
typical (dashed 
line) DOS for a 50 3 
lattice with PBC. 
K 3 x K r = 32 x 32, 
M = 8192. 
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M£)=exp(-^£|iMM£))) • (8) 

On the one hand, the arithmetic mean for large enough K r and K a converges to 
the standard density of states p(E) = Y2n=i ~ E n ), which is not critical at the 
Anderson transition. The geometric mean, on the other hand, represents the typical 
value of the distribution, which, as shown above, is finite in the delocalised phase, 
but goes to zero in the localised phase. As can be seen from Fig. 5, p^{E) and p ty (E) 
are almost equal for extended states, whereas for localised states pt y (E) vanishes 
and pa, v (E) remains finite. This implies, that the ratio of these two quantities, the 
normalised typical density of states 

R(E) := ^ , (9) 

Pav (,-£/) 

can serve as an order parameter for the Anderson transition with R > for extended 
states and R = for localised ones. As for most numerical calculations the transition 
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Fig. 6 Normalised 
typical DOS as a 
function of disorder 
calculated with 
increasing expan- 
sion order on a 50 3 
lattice. The inset 
shows the corre- 
sponding behaviour 
of the Id system 
with 125000 sites. 
K s x K r = 32 x 32. 
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is slightly washed out by the finite size of the considered cluster, and by the KPM 
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resolution. However, for large clusters and increasing M a plot of the ratio 7? versus 
disorder strength W (see Fig. 6) allows for a reliable determination of the critical 
disorder W c , and, e.g., in the band centre we obtain W C {E = 0) ~ 16.54 in accordance 
with other numerical results for the 3d Anderson model [5,32,33]. The quality of 
this criterion is underlined also by our data for a Id system shown in the inset of 
Fig. 6. Here, as mentioned above, arbitrarily small disorder leads to localisation of 
the entire spectrum. Clearly, in our approach this is reflected by a typical DOS that 
vanishes for large M. 

2.4 Comparison of the different methods 

Comparing the value of the critical disorder obtained by the various methods dis- 
cussed in the previous sections, the two main results are the following: (i) As can 
be seen from Fig. 7, the established criteria and methods show an uncertainty of 
the critical value W c in the order of ±0.5t, which is mainly due to the finite system 
sizes accessible to the numerical calculations. Note that our data widely agrees with 
the results in the literature [5,34]. (ii) The value W c — 16.5t can be reproduced 
with the same accuracy using a vanishing typical DOS as an indicator for localisa- 
tion. An improvement of the accuracy of this result can in principle be obtained by 
extending the numerical effort (larger systems, higher resolution, high-performance 
computers) , which is facilitated by the straightforward parallelisability of the KPM 
algorithm. On the other hand, an appropriate scaling ansatz may improve the esti- 
mate of W c on the basis of the presented data. 




Fig. 7 Comparison of the critical values of disorder, W c , obtained by the methods 
outlined in Sees. 2.1-2.3. Decay of the wavefunction (triangles): TV = 20 3 and 30 3 , 
E n G [—0.01,0.01], K r = 10. Average inverse participation number (diamonds): 
N = 16 3 , E n G [-0.1,0.1], K r = 100. Normalised typical DOS (circles): N = 50 3 , 
M = 16384, K r x K a = 32 x 32. 
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Using the well-established value W C (E = 0) ~ 16. 5t as a calibration of the critical 
R, required to distinguish localised from extended states for the used values of N 
and M, we reproduce the mobility edge in the energy-disorder plane [5, 33] using 
R c ~ 0.05 (see Fig. 8). We also find the well-known reentrant behaviour near the 
unperturbed band edges [32,35]: Varying W for some fixed values of E (6i < E < 
7.6t) a region of extended states separates two regions of localised states. The Lifshitz 
boundaries, shown as dashed lines, indicate the energy range, where eigenstates 
are in principle allowed. As the probability of reaching the Lifshitz boundaries is 
exponentially small, we cannot expect to find states near these boundaries for the 
finite ensembles considered in any numerical calculation. 

With respect to numerical resources, clearly, the methods that are based on a 
complete diagonalisation of the Hamiltonian (decay of the wavefunction or partic- 
ipation number) are the least favourable ones, since their CPU requirements scale 
as iV 3 and the memory as iV 2 . Using LAPACK [36] routines for dense matrices on 
a standard PC-system diagonalisations are feasible for systems up to 21 3 . Banded 
matrix routines together with the bandwidth reduction described in Appendix B 
allow to increase this size to 30 3 . 

The calculation of the LDOS via KPM is based on sparse matrix vector multi- 
plications, whose CPU and memory requirements scale only linearly in N. Hence, 
systems up to 100 3 can be easily handled with desktop computers, and the use 
of high-performance environments permits the study of even larger ensembles and 
systems. We conclude that the new method substantially increases the size of nu- 
merically accessible systems, which may lead to a more thorough understanding of 
the Anderson transition. 
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Fig. 8 Phase diagram of the An- 
derson model on a 3d cubic lat- 
tice. Shown are the mobility edge 
(solid curve) as well as the Lif- 
shitz boundaries (dashed lines). 
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3 Conclusions 

With this contribution we aimed to compare well-established numerical localisation 
criteria for electrons in disordered systems with a recently proposed approach that 
is based on the evaluation of the typical density of states. 

Considering the 3d cubic Anderson model we proved that the local DOS can be 
very efficiently calculated using a Chebyshev expansion with kernel polynomial re- 
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finement. Given the numerically obtained distribution of the LDOS, the correspond- 
ing typical DOS allows for the detection of the delocalisation-localisation transition 
with a precision that is comparable to results known from other methods. Like for 
all numerical schemes, the method is restricted to finite systems, and the obtained 
critical values are subjected to finite-size effects. However, the low computational 
resources required by our approach substantially increases the accessible system 
sizes. 

Finally we established the use of the typical DOS as a kind of order parameter, 
which is important in view of its application to interacting disordered systems. 
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A Calculation of the LDOS via the kernel polynomial 
method 

At first glance, Eq. (5) suggests that the calculation of the LDOS could require a 
complete diagonalisation of H. It turns out, however, that an expansion of pi in 
terms of Chebyshev polynomials T n (x) = cos(n arccos x) allows for an incredibly 
precise approximation. Since the Chebyshev polynomials form an orthogonal set on 
the interval [—1, 1], prior to an expansion the Hamiltonian H needs to be rescaled, 

W/2 + 2dt + O.Olt ' V ' 

Here W/2 + 2dt reflects half the bandwidth of the Anderson model and O.Olt is an 
additional offset that ensures numerical stability of the expansion. In terms of the 
coefficients 

r 

Hm= I Pl (x)T m (x)dx = Y^{i\n){n\T m (x n )\i) = (i\T m (X)\i) (11) 

jfj 71 = 1 

the approximate LDOS p~i(x) reads 

1 ( M \ 

pi{x) = — goUo + 2 V g m p m T m (x) . (12) 

The factors 

9m = — ((M - m + 1) cos(m0) + j , (13) 

where = ir/(M + 1), result from a convolution of the finite series with the so- 
called Jackson kernel [17], which mainly damps out the Gibbs oscillations known 
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from polynomial approximations (cf. Fig. 9). The width of the kernel, Ax = tt/M, 
scales inversely with the order of the expansion M and defines the resolution of the 
method. 
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Fig. 9 Chebyshev expansion of a 8- 
peak: The plain truncated series of 
order M = 64 is a strongly oscillating 
curve (dashed). By convolution with 
the Jackson kernel it transforms into 
a strictly positive, well localised peak 
at x — (solid) , which is much closer 
to our usual notion of S(x). 



Using the recursion relations of the Chebyshev polynomials, 

T m+1 (x) = 2xT m {x) - T m -i(x) , (14) 

the moments [i m can be calculated iteratively. An additional trick allows for the 
generation of two moments with each matrix vector multiplication by X, 

N 

fJtom-i = ^2 2{i\T m (X)T m - 1 (X)\i) - fix , 

1=1 (15) 

M2m = 2(i|T m (X)r m (X)jj) — fio , 
i=l 

reducing the numerical effort by another factor 1/2. Note that the algorithm requires 
storage only for the sparse matrix X and two vectors of the corresponding dimension. 



B Reduction of the bandwidth of the Anderson matrix 
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Fig. 10 Sparsity pattern of the An- 
derson matrix. Left: Standard An- 
derson matrix for a 9 3 system with 
PBC. The use of DSBEV would 
require the storage of 648 off- 
diagonals. Right: Reduced matrix af- 
ter transformation (16). Only 162 off- 
diagonals have to be stored. 



While for Krylov sub-space methods [37] (like Lanczos or Jacobi-Davidson) we 
can take advantage of the sparsity of the tight-binding type matrices, the full diag- 
onalisation with LAPACK routines requires their complete storage. Unfortunately, 
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for periodic boundary conditions the cyclic tridiagonal structure of the matrices 
spoils the use of band matrix routines like DSBEV. For a L 3 -system there are non- 
vanishing matrix elements in a distance of L z — L 2 from the diagonal (Fig. 10). 
Thus almost all matrix elements (most of them zero) need to be stored, giving no 
advantage compared to full matrix diagonalisation routines like DSYEV. 

For linear systems there are tricks to use tridiagonal matrices instead of the cor- 
responding cyclic tridiagonal matrices, which are based on the use of the Sherman- 
Morrison formula [38] . We are not aware of similar ideas for eigenvalue problems. It 
turns out, however, that an appropriate sequence of Givens rotations [38] allows the 
transformation of the cyclic tight-binding matrix (with quadratic blocks close to the 
outer edges) onto a matrix with blocks only along five diagonals. The corresponding 
transformation reads 

H Tcd = T T HT, (16) 
where T = P ® l L 2 xL 2 and for odd L the L x L matrix P is given by 



P = — 




(17) 



For even L the first row and column are absent. 

The bandwidth of H can thus be substantially reduced (to 2L 2 , see Fig. 10), 
which allows for the full diagonalisation of systems up to 30 3 sites on PC-systems 
with a memory of 512 MB. Furthermore, the sparsity of the transformation (16) can 
be used to avoid an explicit matrix-matrix multiplication [39]. Hence the change 
from H to H le d is not time consuming. The advantage of the transformation is not 
primarily a gain in CPU time but storage. 
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